Function Optimization

1 Unconstrained optimization

In this document, we explore a series of optimization problems ranging from simple unconstrained optimization to more complex denoising tasks. We employ a variety of R packages and techniques to solve and visualize these problems. ## Problem 1

# Define the function
f1 <- function(x) {
  return(x^4)
}

# Apply optim() to find the minimum
result1 <- optim(0, f1, method = "BFGS") # Starting from 0 as an initial guess
result1$par # Optimal value of x

## [1] 0

result1$value # Minimum value of f(x)

## [1] 0

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.3

library(plotly)

## Warning: package 'plotly' was built under R version 4.2.3

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# Plot the function
x_vals <- seq(-2, 2, by=0.01)
y_vals <- x_vals^4
df1 <- data.frame(x=x_vals, y=y_vals)

ggplot(df1, aes(x=x, y=y)) +
  geom_line(color="blue") +
  ggtitle(expression(f(x) == x^4)) +
  theme_minimal()

The minimum value is achieved at x=0. The minimum value of the function is 0.

1.1 Problem 2

# Define the function
f2 <- function(x) {
  return(-(2*sin(x) - x^2/10)) # Negate the function for maximization
}

# Apply optim() to find the maximum
result2 <- optim(0, f2, method = "BFGS") # Starting from 0 as an initial guess
result2$par # Optimal value of x

## [1] 1.427552

result2$value # Maximum value of f(x) (Note: this will be the negated value)

## [1] -1.775726

# Plot the function
y_vals <- 2*sin(x_vals) - x_vals^2/10
df2 <- data.frame(x=x_vals, y=y_vals)

ggplot(df2, aes(x=x, y=y)) +
  geom_line(color="red") +
  ggtitle(expression(f(x) == 2*sin(x) - x^2/10)) +
  theme_minimal()

The maximum value is achieved at x≈1.427552. The maximum value of the function is approximately −1.775726.

1.2 Problem 3

# Define the function
f3 <- function(params) {
  x <- params[1]
  y <- params[2]
  return(-(2*x*y + 2*x - x^2 - 2*y^2)) # Negate the function for maximization
}

# Apply optim() to find the maximum
initial_guess <- c(0,0) # Initial guess for x and y
result3 <- optim(initial_guess, f3, method = "BFGS") 
result3$par # Optimal values of x and y

## [1] 2.000003 1.000002

result3$value # Maximum value of f(x, y) (Note: this will be the negated value)

## [1] -2

# Create a 3D plot

x_vals <- seq(-10, 10, by=0.5)
y_vals <- seq(-10, 10, by=0.5)
z <- matrix(0, length(x_vals), length(y_vals))

for (i in 1:length(x_vals)) {
  for (j in 1:length(y_vals)) {
    z[i,j] <- 2*x_vals[i]*y_vals[j] + 2*x_vals[i] - x_vals[i]^2 - 2*y_vals[j]^2
  }
}

persp(x_vals, y_vals, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue", 
      xlab = "x", ylab = "y", zlab = "f(x,y)", main = expression(f(x, y) == 2*x*y + 2*x - x^2 - 2*y^2))

The maximum value is achieved at the point (x,y)≈(2.000003,1.000002). The maximum value of the function is −2.

2 Linear Programming

2.1 Using lpSolveAPI

# Define the objective function coefficients
f.obj <- c(1, 2, 3, 4)

# Define the constraint matrix
f.con <- matrix(c(4, 3, 2, 1,
                 1, 0, -1, 2,
                 1, 1, 1, 1), 
                nrow=3, byrow=TRUE)

# Define the direction of the constraints
f.dir <- c("=", ">=", ">=")

# Define the right-hand side of the constraints
f.rhs <- c(2, 1, 0)

# Create a linear program model
lp.model <- make.lp(0, 4)

# Set the objective function coefficients
set.objfn(lp.model, f.obj)

# Add constraints
add.constraint(lp.model, f.con[1,], f.dir[1], f.rhs[1])
add.constraint(lp.model, f.con[2,], f.dir[2], f.rhs[2])
add.constraint(lp.model, f.con[3,], f.dir[3], f.rhs[3])

# Set bounds for decision variables
set.bounds(lp.model, column = 1, lower = 0)
set.bounds(lp.model, column = 3, lower = 0)
set.bounds(lp.model, column = 4, upper = 10)

# Set the direction of optimization (maximization)
lp.control(lp.model, sense="max")

## $anti.degen
## [1] "fixedvars" "stalling" 
## 
## $basis.crash
## [1] "none"
## 
## $bb.depthlimit
## [1] -50
## 
## $bb.floorfirst
## [1] "automatic"
## 
## $bb.rule
## [1] "pseudononint" "greedy"       "dynamic"      "rcostfixing" 
## 
## $break.at.first
## [1] FALSE
## 
## $break.at.value
## [1] 1e+30
## 
## $epsilon
##       epsb       epsd      epsel     epsint epsperturb   epspivot 
##      1e-10      1e-09      1e-12      1e-07      1e-05      2e-07 
## 
## $improve
## [1] "dualfeas" "thetagap"
## 
## $infinite
## [1] 1e+30
## 
## $maxpivot
## [1] 250
## 
## $mip.gap
## absolute relative 
##    1e-11    1e-11 
## 
## $negrange
## [1] -1e+06
## 
## $obj.in.basis
## [1] TRUE
## 
## $pivoting
## [1] "devex"    "adaptive"
## 
## $presolve
## [1] "none"
## 
## $scalelimit
## [1] 5
## 
## $scaling
## [1] "geometric"   "equilibrate" "integers"   
## 
## $sense
## [1] "maximize"
## 
## $simplextype
## [1] "dual"   "primal"
## 
## $timeout
## [1] 0
## 
## $verbose
## [1] "neutral"

# Solve the LP
solve(lp.model)

## [1] 0

# Get the solution
get.variables(lp.model) # check against the exact solution x_1 = 0, x_2 = 0, x_3 = 0, x_4=2

## [1] 0 0 0 2

get.objective(lp.model) #get optimal max value

## [1] 8

# Define the constraints
x1 <- seq(0, 10, length.out = 100)
df <- data.frame(x1 = x1)

# Using the first constraint 4x1 + 3x2 + 2x3 + x4 = 2 and fixing x3=0, x4=0
df$x2_constraint1 <- (2 - 4*df$x1)/3

# Using the second constraint x1 - x3 + 2x4 >= 1 and fixing x3=0, x4=0
df$x2_constraint2 <- Inf  # This constraint doesn't give a relation between x1 and x2 for x3=0, x4=0

# Using the third constraint x1 + x2 + x3 + x4 >= 0 and fixing x3=0, x4=0
df$x2_constraint3 <- -df$x1

# Plotting the feasible region
ggplot(df, aes(x1)) +
  geom_line(aes(y = x2_constraint1), color = "blue", linetype = "dashed") +
  geom_line(aes(y = x2_constraint2), color = "red", linetype = "dashed") +
  geom_line(aes(y = x2_constraint3), color = "green", linetype = "dashed") +
  # Highlighting the feasible region
  geom_ribbon(aes(ymin = pmax(x2_constraint1, x2_constraint3), ymax = 10), fill = "grey80") +
  labs(title = "Feasible Region in x1-x2 Plane",
       x = "x1", y = "x2") +
  theme_minimal()

The solve(lp.model) result of 0 indicates that a feasible solution was found and the optimization was successful. In lpSolveAPI, a return value of 0 means the problem has an optimal solution.

get.variables(lp.model) returns the optimal values of the decision variables in the LP problem. The result 0 0 0 2 corresponds to x1=0, x2=0, x3=0, and x4=2. This means that to maximize the objective function given the constraints, you should set x4 to 2 and keep the other variables at 0.

get.objective(lp.model) returns the optimal value of the objective function given the optimal solution. The result 8 indicates that the maximum value of the objective function (given the constraints) is 8.

2.2 Using Rsolnp

# Objective function
objective_func <- function(x) {
  return(-1 * (x[1] + 2*x[2] + 3*x[3] + 4*x[4] + 5))  # We multiply by -1 because solnp minimizes by default
}

# Constraints
constraints <- function(x) {
  c(
    4*x[1] + 3*x[2] + 2*x[3] + x[4] - 2,
    x[1] - x[3] + 2*x[4] - 1,
    x[1] + x[2] + x[3] + x[4]
  )
}

# Initial values
x0 <- c(0, 0, 0, 0)

# Solve using Rsolnp
sol <- solnp(
  pars = x0,
  fun = objective_func,
  eqfun = constraints,
  eqB = c(0, 0, 0),
  LB = c(0, -Inf, 0, -Inf),
  UB = c(Inf, Inf, Inf, 10)
)

## Warning in p0 * vscale[(neq + 2):(nc + np + 1)]: longer object length is not a
## multiple of shorter object length

## 
## Iter: 1 fn: -5.0000   Pars:  0       0       0       0      
## solnp--> Solution not reliable....Problem Inverting Hessian.

solution_Rsolnp <- sol$pars
solution_Rsolnp

## [1] 0 0 0 0

The solution is not reliable due to a problem inverting the Hessian.The repeated warning solnp means the linearized problem has no feasible solution. Solution for the variables x1,x2,x3,x4: [0, 0, 0, 0]

Comparison: While lpSolveAPI found a feasible and optimal solution, solnp struggled and returned infeasible values for the decision variables.

lpSolveAPI found an optimal objective value of 8, while solnp didn’t provide a clear optimal value due to its issues.

lpSolveAPI suggests x4=2 (and other variables at 0) as the optimal decision, while solnp provides unrealistic values.

In Conclusion: For this particular LP problem, lpSolveAPI provided a feasible and optimal solution, while solnp faced challenges. It’s important to note that while solnp is a versatile optimizer capable of handling non-linear problems, it might not be the best choice for linear problems, especially when there are specific solvers like lpSolveAPI designed for LPs.

3 Mixed Integer Linear Programming

# Create a new LP model with 2 decision variables
lprec <- make.lp(0, 2)

# Set the objective function
# min: 4*x1 + 6*x2
set.objfn(lprec, c(4, 6))

# Add constraints
# 2*x1 + 2*x2 >= 5
add.constraint(lprec, c(2, 2), ">=", 5)
# x1 - x2 <= 1
add.constraint(lprec, c(1, -1), "<=", 1)

# Define bounds on the decision variables
# Both x1 and x2 are >= 0
set.bounds(lprec, lower = rep(0, 2), upper = rep(Inf, 2))

# Set x1 and x2 to be integer variables
set.type(lprec, columns = c(1,2), type = "integer")

# Set direction to minimize
lp.control(lprec, sense = "min")

## $anti.degen
## [1] "fixedvars" "stalling" 
## 
## $basis.crash
## [1] "none"
## 
## $bb.depthlimit
## [1] -50
## 
## $bb.floorfirst
## [1] "automatic"
## 
## $bb.rule
## [1] "pseudononint" "greedy"       "dynamic"      "rcostfixing" 
## 
## $break.at.first
## [1] FALSE
## 
## $break.at.value
## [1] -1e+30
## 
## $epsilon
##       epsb       epsd      epsel     epsint epsperturb   epspivot 
##      1e-10      1e-09      1e-12      1e-07      1e-05      2e-07 
## 
## $improve
## [1] "dualfeas" "thetagap"
## 
## $infinite
## [1] 1e+30
## 
## $maxpivot
## [1] 250
## 
## $mip.gap
## absolute relative 
##    1e-11    1e-11 
## 
## $negrange
## [1] -1e+06
## 
## $obj.in.basis
## [1] TRUE
## 
## $pivoting
## [1] "devex"    "adaptive"
## 
## $presolve
## [1] "none"
## 
## $scalelimit
## [1] 5
## 
## $scaling
## [1] "geometric"   "equilibrate" "integers"   
## 
## $sense
## [1] "minimize"
## 
## $simplextype
## [1] "dual"   "primal"
## 
## $timeout
## [1] 0
## 
## $verbose
## [1] "neutral"

# Solve the LP
solve(lprec)

## [1] 0

# Extract solution
solution <- get.variables(lprec)

# Print the solution
solution

## [1] 2 1

# Constraints
x1 <- seq(0, 5, by = 0.1)
x2_1 <- (5 - 2*x1) / 2  # 2x1 + 2x2 >= 5
x2_2 <- x1 - 1          # x1 - x2 <= 1

# Plot
ggplot(data.frame(x1 = x1), aes(x1)) +
  geom_line(aes(y = x2_1), color = "blue", linetype="dashed") +
  geom_line(aes(y = x2_2), color = "red", linetype="dashed") +
  annotate("text", x = 3, y = 3, label = "2x1 + 2x2 >= 5", color="blue") +
  annotate("text", x = 3, y = 1, label = "x1 - x2 <= 1", color="red") +
  geom_polygon(aes(y = pmin(x2_1, x2_2)), fill = "gray80") + 
  geom_point(aes(x = 2, y = 1), color = "green", size = 4) +  # Optimal solution
  labs(title = "Feasible Region and Optimal Solution for MILP",
       x = expression(x[1]),
       y = expression(x[2])) +
  theme_minimal()

The output indicates that the solver successfully found a solution (a return value of 0 means the problem was successfully solved). The optimal solution is x1=2,x2=1 for the objective function 4x1+6x^2 to be minimized.

4 Quadratic Programming

4.1 Using quadprog

library(quadprog)
# Define the matrices and vectors for the QP
Dmat <- matrix(c(4, 1, 
                 1, 2), 2, 2) # Coefficient matrix for the quadratic terms
dvec <- c(-1, -1) # Coefficient vector for the linear terms

# Coefficient matrix for the constraints
Amat <- matrix(c(-1, 1, 0,
                 -1, 0, 1), ncol=2) 
bvec <- c(-1, 0, 0) # Right-hand side of the constraints

# Solve the QP
solution_qp <- solve.QP(Dmat, dvec, t(Amat), bvec)

# Extract the solution
x_qp <- solution_qp$solution
solution_qp$value

## [1] 2.775558e-17

Given the quadratic objective function and constraints provided, quadprog determined that the minimum of the objective function occurs when both x1 and x2 are 0. This means, within the constraints provided, the function reaches its minimum value when both variables are set to zero

4.2 Using Rsolnp

library(Rsolnp)

# Objective function
obj_fun <- function(x) {
  return(2*x[1]^2 + x[2]^2 + x[1]*x[2] + x[1] + x[2])
}

# Equality constraint
eq_fun <- function(x) {
  return(x[1] + x[2] - 1)
}

# Initial guess
x0 <- c(0.5, 0.5)

# Solve using Rsolnp
sol <- solnp(pars = x0, fun = obj_fun, eqfun = eq_fun, eqB = 0)

## 
## Iter: 1 fn: 1.8750    Pars:  0.25000 0.75000
## Iter: 2 fn: 1.8750    Pars:  0.25000 0.75000
## solnp--> Completed in 2 iterations

# Extract solution
sol$pars

## [1] 0.25 0.75

sol$values[length(sol$values)]

## [1] 1.875

Given the same objective function and constraints, Rsolnp found that the function reaches its minimum when x1 is 0.25 and x2 is 0.75.

This suggests that there might be some differences in how the two packages handle the optimization problem, or it could be due to the nonlinear nature of the solver introducing some approximation.

4.3 Determine the Automated Lagrange Multiplier(Rsolnp)

The Automated Lagrange multiplier can be extracted directly from the Rsolnp result.

lagrange_qp <- sol$lagrange
lagrange_qp

##      [,1]
## [1,] 2.75

The Lagrange multiplier provides insights about the constraints. A positive value of the Lagrange multiplier (like 2.75) suggests that if we increase the right-hand side of the equality constraint by a small amount, the objective function would increase by the value of the Lagrange multiplier (assuming a maximization problem).

In simpler terms, the constraint is binding, and the solution is sensitive to changes in the constraint’s limit.

4.4 Manual Lagrange multiplier

To find the Lagrange multipliers manually, we’ll set up the Lagrangian and then find its gradient with respect to both the original variables and the multipliers. Setting this gradient to zero will give us the necessary equations.

library(numDeriv)
# Lagrangian function with adjusted sign for constraint
lagrangian_func<- function(params) {
  x1 <- params[1]
  x2 <- params[2]
  lambda <- params[3]
  
  # Objective function
  obj <- 2*x1^2 + x2^2 + x1*x2 + x1 + x2
  
  # Constraint with adjusted sign
  constraint <- -(x1 + x2 - 1)
  
  # Lagrangian with adjusted constraint
  L <- obj + lambda * constraint
  return(L)
}

# Calculate gradient and Hessian for the corrected Lagrangian
grad_and_hess <- hessian(lagrangian_func, c(0.5, 0.5, 0))

# Solve for the Lagrange multiplier manually using the corrected Lagrangian
sol_numDeriv <- solve(grad_and_hess, -grad(lagrangian_func, c(0.5, 0.5, 0)))

# Extracting the Lagrange multiplier from the solution
lagrange_numDeriv<- sol_numDeriv[3]

lagrange_numDeriv

## [1] 2.75

The manually computed Lagrange multiplier matches the automated result from Rsolnp, which is a good validation. This value further confirms the importance and sensitivity of the constraint to the solution of the problem.

4.5 Compare the three versions of the results

To find the Lagrange multipliers manually, we’ll set up the Lagrangian and then find its gradient with respect to both the original variables and the multipliers.

Setting this gradient to zero will give us the necessary equations.

cat("Solution using quadprog:", x_qp, "\n")

## Solution using quadprog: 0 0

cat("Solution using Rsolnp:", sol$pars, "\n")

## Solution using Rsolnp: 0.25 0.75

cat("Automated Lagrange Multiplier(Rsolnp):", lagrange_qp, "\n")

## Automated Lagrange Multiplier(Rsolnp): 2.75

cat("Manual Lagrange multiplier", lagrange_numDeriv, "\n")

## Manual Lagrange multiplier 2.75

The solutions from quadprog and Rsolnp differ. This could be due to differences in the algorithms or methods used by each package or the way the problem is represented internally.

The Lagrange multipliers, both automated and manual, match and suggest that the constraint is binding. This means that the solution is on the boundary of the feasible region defined by the constraint, and any change in the constraint will affect the optimal solution.

# Load necessary libraries
library(ggplot2)

# Define the constraints
x1 <- seq(0, 1, length.out = 300)
x2 <- 1 - x1

# Define the objective function for coloring
obj_func <- function(x1, x2) {
  return(2*x1^2 + x2^2 + x1*x2 + x1 + x2)
}

# Create a data frame for plotting
df <- data.frame(x1=x1, x2=x2, obj_val=obj_func(x1, x2))

# Optimal solution from sol$pars
opt_x1 <- 0.25
opt_x2 <- 0.75


# Define a grid of x1 and x2 values
x1_vals <- seq(0, 1.5, length.out = 300)
x2_vals <- seq(0, 1.5, length.out = 300)
grid <- expand.grid(x1 = x1_vals, x2 = x2_vals)
grid$obj_val <- with(grid, 2*x1^2 + x2^2 + x1*x2 + x1 + x2)

# Plot
ggplot(grid, aes(x=x1, y=x2)) +
  geom_tile(aes(fill=obj_val)) +
  geom_contour(aes(z=obj_val), color="white") +
  geom_point(aes(x=opt_x1, y=opt_x2), color="red", size=4) +
  scale_fill_gradient2(low="blue", mid="white", high="red", midpoint=median(grid$obj_val)) +
  labs(title="Objective Function Contour Plot",
       fill="Objective\nFunction Value") +
  theme_minimal()

The feasible region is visualized with the objective function value varying from blue (low values) to red (high values). The optimal solution is marked with a red point at coordinates (0.25, 0.75).

5 Complex non-linear optimization

# Define the Rosenbrock function
rosenbrock <- function(x) {
  100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
}

# Initial guess
x0 <- c(0.5, 0.5)

# Define the constraints
# Constraint function for x1, x2 >= 0
constraint_fun <- function(x) {
  c(x[1], x[2])
}
constraint_lb <- c(0, 0)
constraint_ub <- c(Inf, Inf)  # Since there are no upper bounds on x1 and x2

# Solve the optimization problem
sol <- solnp(pars = x0, fun = rosenbrock, ineqfun = constraint_fun, ineqLB = constraint_lb, ineqUB = constraint_ub)

## 
## Iter: 1 fn: 9.29e-10  Pars:  0.99997 0.99994
## Iter: 2 fn: 9.035e-10     Pars:  0.99997 0.99994
## solnp--> Completed in 2 iterations

# Display the solution
sol$values

## [1] 6.500000e+00 9.289547e-10 9.034685e-10

sol$pars

## [1] 0.9999699 0.9999398

# Load necessary libraries
library(lattice)

# Define the function f(x1, x2)
f <- function(x1, x2) {
  return(100*(x2 - x1^2)^2 + (1 - x1)^2)
}

# Create a grid of x1 and x2 values
x1_seq <- seq(0, 2, length.out = 100)
x2_seq <- seq(0, 2, length.out = 100)
grid <- expand.grid(x1 = x1_seq, x2 = x2_seq)

# Calculate f(x1, x2) for each grid point
grid$z <- with(grid, f(x1, x2))

# Plot the surface
wireframe(z ~ x1 * x2, data = grid, drape = TRUE, colorkey = TRUE, 
          scales = list(arrows = FALSE),
          zlab = "f(x1, x2)", xlab = "x1", ylab = "x2",
          main = "Function Surface with Optimal Solution")

The sol$values are very close to zero, indicating that the algorithm is converging to the minimum value of the function.

In this run, the solnp algorithm started from an initial guess of x0 = c(0.5, 0.5) and within 2 iterations, it has converged very close to the known minimum at x1 = 1 and x2 = 1. The values obtained are x1 = 0.9999699 and x2 = 0.9999398, which are extremely close to 1, confirming the efficacy of the algorithm for this problem.

In essence, the algorithm did an excellent job of finding the known global minimum of the Rosenbrock function (1,1) subject to the given constraints.

6 Data Denoising

denoise_data <- function(noise_level) {
  library(plotly)
  library(Rsolnp)
  library(devtools)
  #install_github("cran/tvd")  # 
  library("tvd")
  
  n <- 1000
  xs <- seq(0, 8, len=n)
  
  x_noisy <- function (x) { 
    sin(x)^2/(1.5+cos(x)) + rnorm(length(x), 0, noise_level)
  }


set.seed(1234)

x_noisy <- function (x) { 
  sin(x)^2/(1.5+cos(x)) + rnorm(length(x), 0, noise_level)
}

# initialize the manual denoised signal
x_denoisedManu <- rep(0, n)

df <- as.data.frame(cbind(xs, x_noisy(xs)))
# loess fit a polynomial surface determined by numerical predictors,
# using local fitting
poly_model1 <- loess(x_noisy(xs) ~ xs, span=0.1, data=df) # tight model
poly_model2 <- loess(x_noisy(xs) ~ xs, span=0.9, data=df) # smoother model
# To see some of numerical results of the model-fitting:
# View(as.data.frame(cbind(xs, x_noisy, predict(poly_model1))))

# plot(xs, x_noisy(xs), type='l')
# lines(xs, poly_model1$fitted, col="red", lwd=2)
# lines(xs, poly_model2$fitted, col="blue", lwd=3)

plot_ly() %>%
  add_trace(x=~xs, y=~x_noisy(xs), type="scatter", mode="lines", name="Noisy Data") %>%
  add_trace(x=~xs, y=~poly_model1$fitted, type="scatter", mode="lines", line=list(width=5), name="Tight LOESS Smoothing") %>%
  add_trace(x=~xs, y=~poly_model2$fitted, type="scatter", mode="lines", line=list(width=5), name="Smoother LOESS Model") %>%
  add_trace(x=~xs, y=~sin(xs)^2/(1.5+cos(xs)), type="scatter", mode="lines", line=list(width=5), name="Original Process") %>%
  layout(title="Noisy Data and Smoothed Models", legend = list(orientation='h'))
#Next, let’s initiate the parameters, define the objective function and optimize it, i.e., estimate the parameters that minimize the cost function as a mixture of fidelity and regularization terms.

# initialization of parameters
betas_0 <- c(0.3, 0.3, 0.5, 1)
betas <- betas_0

# Denoised model
x_denoised <- function(x, betas) {
  if (length(betas) != 4) {
    print(paste0("Error!!! length(betas)=", length(betas), " != 4!!! Exiting ..."))
    break();
  }  
  # print(paste0("   .... betas = ", betas, "...\n"))
  # original noise function definition: sin(x)^2/(1.5+cos(x))
  return((betas[1]*sin(betas[2]*x)^2)/(betas[3]+cos(x)))
} 

library(Rsolnp)

fidelity <- function(x, y) {
  sqrt((1/length(x)) * sum((y - x)^2))
}

regularizer <- function(betas) {
  reg <- 0
  for (i in 1:(length(betas-1))) {
      reg <- reg + abs(betas[i])
  }
  return(reg)
}

# Objective Function
objective_func <- function(betas)  {
# f(x) = 1/2 * \sum_{t=1}^{n-1} {|y(t) - x_{noisy}(t)\|^2}} + \lambda * \sum_{t=2}^{n-1} | x(t) - x(t-1)|
  fid <- fidelity(x_noisy(xs), x_denoised(xs, betas)) 
  reg <- abs(betas[4])*regularizer(betas)
  error <- fid + reg
  # uncomment to track the iterative optimization state
  # print(paste0(".... Fidelity =", fid, " ... Regularizer = ", reg, " ... TotalError=", error))
  # print(paste0("....betas=(",round(betas[1],3),",",round(betas[2],3),",",
  #       round(betas[3],3),",", round(betas[4],3), ")"))
  return(error)
}

# inequality constraint forcing the regularization parameter lambda=beta[4]>0
inequalConstr <- function(betas){
   betas[4]
}
inequalLowerBound <- 0; inequalUpperBound <- 100

# should we list the value of the objective function and the parameters at every iteration (default trace=1; trace=0 means no interim reports)
# constraint problem 
ctrl <- list(trace=0, tol=1e-5)  ## could specify: outer.iter=5, inner.iter=9)
set.seed(121)
sol_lambda <- solnp(betas_0, fun = objective_func, ineqfun = inequalConstr, ineqLB = inequalLowerBound, ineqUB = inequalUpperBound, control=ctrl)

# unconstrained optimization
# ctrl <- list(trace=1, tol=1e-5)  ## could specify: outer.iter=5, inner.iter=9)
# sol_lambda <- solnp(betas_0, fun = denoising_func, control=ctrl)

# suppress the report of the functional values (too many)
# sol_lambda$values

# report the optimal parameter estimates (betas)
sol_lambda$pars
## [1] 2.5649689 0.9829681 1.7605481 0.9895268
# Reconstruct the manually-denoised signal using the optimal betas
betas <- sol_lambda$pars
x_denoisedManu <- x_denoised(xs, betas)
print(paste0("Final Denoised Model:", round(betas[1],3), "*sin(", 
             round(betas[2],3), "*x)^2/(", round(betas[3],3), "+cos(x)))"))
## [1] "Final Denoised Model:2.565*sin(0.983*x)^2/(1.761+cos(x)))"
# plot(x_denoisedManu)

plot_ly() %>%
  add_trace(x=~xs, y=~x_noisy(xs), type="scatter", mode="lines", name="Noisy Data") %>%
  add_trace(x=~xs, y=~x_denoisedManu/2, type="scatter", mode="lines", name="Manual denoising", line=list(width=5)) %>%
  add_trace(x=~xs, y=~sin(xs)^2/(1.5+cos(xs)), type="scatter", mode="lines", line=list(width=5), name="Original Process") %>%
  layout(title="Noisy Data and Smoothed Models", legend = list(orientation='h'))
#Finally, we can validate our manual denoising protocol against the automated TV denoising using the R package tvd, may need to install the RTools package which allows building r-packages from source code (including the linux make tool).


lambda_0 <- 0.5
x_denoisedTVD <- tvd1d(x_noisy(xs), lambda_0, method = "Condat")

# lambda_o is the total variation penalty coefficient
# method is a string indicating the algorithm to use for denoising. 
# Default method is "Condat"


 # Return the plotly object for further use or rendering
  p <- plot_ly() %>%
    add_trace(x=~xs, y=~x_noisy(xs), type="scatter", mode="lines", name="Noisy Data", opacity=0.3) %>%
    add_trace(x=~xs, y=~x_denoisedManu/2, type="scatter", mode="lines", name="Manual Denoising", line=list(width=5)) %>%
    add_trace(x=~xs, y=~poly_model1$fitted, type="scatter", mode="lines", name="LOESS Smoothing", line=list(width=5)) %>%
    add_trace(x=~xs, y=~sin(xs)^2/(1.5+cos(xs)), type="scatter", mode="lines", line=list(width=5), name="Original Process") %>%
    add_trace(x=~xs, y=~x_denoisedTVD, type="scatter", mode="lines", name="TV Denoising", line=list(width=5)) %>%
    layout(title="Noisy Data and Smoothed Models", legend = list(orientation='h'))
  
  return(p)
}

# Test the function (random increment of noise level)
denoise_data(0.111)

## Loading required package: usethis

## [1] "Final Denoised Model:2.18*sin(-0.867*x)^2/(1.516+cos(x)))"

Noise Level: 0.111 At this relatively low noise level, the denoising model seems to fit the original process quite well, capturing the general trend and maintaining smoothness.The manual denoising model has coefficients that present a relatively straightforward transformation of the sine function.

denoise_data(0.232)

## [1] "Final Denoised Model:2.811*sin(1.01*x)^2/(1.792+cos(x)))"

denoise_data(0.597)

## [1] "Final Denoised Model:2.238*sin(1.012*x)^2/(1.789+cos(x)))"

Noise Level: 0.233 & 0.597 For these noise levels, the manual denoising adjusts itself to fit the original process quite closely. The TV denoising, especially at the 0.597 noise level, becomes a smoother representation.

denoise_data(0.342)

## [1] "Final Denoised Model:0.672*sin(-0.194*x)^2/(1.282+cos(x)))"

Noise Level: 0.342 At this particular noise level, both the manual denoising and the TV denoising exhibit deviations from the original process. The manual denoising, while attempting to fit the original process, has some discrepancies, especially at the peaks. The TV denoising appears to be under-smoothed, and it doesn’t capture the pattern of the original process effectively.

denoise_data(0.786)

## [1] "Final Denoised Model:4.489*sin(-0.104*x)^2/(2.117+cos(x)))"

denoise_data(0.903)

## [1] "Final Denoised Model:4.301*sin(-0.017*x)^2/(3.221+cos(x)))"

Noise Level: 0.786 & 0.903 At these higher noise levels, the challenge of fitting the original process becomes more evident. While the manual denoising still provides a reasonable fit, it begins to deviate slightly from the original process, especially at the peaks and troughs. The TV denoising appears more smoothed out.

Trends based on the Figures: Noise Level Impact: As the noise level increases (from 0.1 to 0.9), the noisy data (in blue) becomes more dispersed and more challenging to fit accurately. LOESS Smoothing: The LOESS Smoothing (in green) tends to fit the original process (in red) quite closely across different noise levels. The model parameters change with varying noise levels to achieve this fit. TV Denoising: The TV (Total Variation) denoising (in purple) also tries to fit the original process but may produce a smoother or sometimes under-smoothed signal depending on the noise level. Manual Denoising: The Manual Denoising (in orange) produces varying fits depending on its span parameter and the noise level. In some cases, it tends to over-smooth.

TV Denoising as Noise Level Changes:

The TV denoising method aims to preserve sharp transitions in data (like edges in images) and remove noise. At lower noise levels (0.1 - 0.3), TV denoising is quite close to the original process. However, as noise levels increase (0.4 and above), TV denoising tends to produce a more smoothed version of the signal.

This smoothing might lead to a loss in some finer details of the original process. At very high noise levels (0.7 - 0.9), TV denoising might not capture sudden changes in the signal as effectively.

Overall Comparison:

Manual denoising provides a flexible method to adapt to different noise levels, offering a close fit to the original process across a range of noise levels.However, as noise levels increase significantly, the fidelity of the fit might reduce slightly.

TV denoising offers a robust method to smooth out noise, especially preserving sharp transitions. However, at high noise levels, it may tend to over-smooth, losing some nuances of the original process.

Through a series of optimization tasks, we’ve seen the versatility of R in tackling both simple and complex problems. The discrepancies in some methods underscore the importance of method selection based on problem characteristics. Future steps might involve exploring other optimization techniques or diving deeper into parameter tuning.